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We construct cosmological long-wavelength solutions without symmetry in general gauge con¬ 
ditions which are compatible with the long-wavelength scheme. We then specify the relationship 
among the solutions in different time slicings. Nonspherical long-wavelength solutions are partic¬ 
ularly important for primordial structure formation in the epoch of very soft equations of state. 

Applying this general framework to spherical symmetry, we show the equivalence between long- 
wavelength solutions in the constant mean curvature slicing with conformally flat spatial coordinates 
and asymptotic quasihomogeneous solutions in the comoving slicing with the comoving threading. 

We derive the correspondence relation between these two solutions and compare the results of nu¬ 
merical simulations of primordial black hole (PBH) formation in these two different approaches. To 
discuss the PBH formation, it is convenient and conventional to use 5c, the value which the averaged 
density perturbation at threshold in the comoving slicing would take at horizon entry in the lowest- 
order long-wavelength expansion. We numerically find that within (approximately) compensated 
models, the sharper the transition from the overdense region to the Friedmann-Robertson-Walker 
universe is, the larger the 5c becomes. We suggest that, for the equation of state p = (P — l)p, we 
can apply the analytic formulas for the minimum 5c,min — [3r /(SP-1-2)] sin^ [ttVP — 1/(3P — 2)] and 
the maximum 5c,max — 3P/(3P + 2). As for the threshold peak value of the curvature variable ipo.c, 
we find that the sharper the transition is, the smaller the ipo^c becomes. We analytically explain this 
intriguing feature qualitatively with a compensated top-hat density model. Plsing simplified models, 
we also analytically deduce an environmental effect that ipo^c can be significantly larger (smaller) if 
the underlying density perturbation of much longer wavelength is positive (negative). 
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I. INTRODUCTION 


The proposal that black holes may have formed in the early Universe, which are called primordial black holes 
(PBHs), has recently been extensively studied in cosmology because they can convey unique information about the 
early Universe to us. This is because energy scales relevant to PBHs are generally very different from those of the 
observed cosmic microwave background anisotropy. The possibility of PBHs has been indicated by Zeldovich and 
Novikov [I| and Hawking 0. Since PBHs evaporate by emitting radiation and they also act as gravitational sources 
before they evaporate away, current observations constrain the abundance of PBHs using the big bangmucleosynthesis, 
extragalactic photon background and gravitational and astrophysical effects of nonevaporating PBHs To convert 

the observed abundance of PBHs to the information of the early Universe, the key issues are the production efficiency 
of PBHs from ^en cosmological perturbations and the time evolution of PBHs through mass accretion and Hawking 
evaporation [11^1!Hi- 

Here, we focus on the criterion of PBH formation. Since the probability of PBH formation is expected to be 
exponentially small, the abundance of PBHs for given perturbations is sensitive to the formation criterion. As a 
pioneer, Carr Q derived i5c — P — 1 for the equation of state p = (P — l)p almost in an order-of-magnitude estimate, 
where 5c is the threshold value of the density perturbation at horizon entry. This gives 5c — 1/3 for a radiation 
fluid r = 4/3. In fact, it turns out that it is not so straightforward to accurately determine the formation criterion 
of PBHs as we have to understand highly general relativistic nonlinear dynamics in the cosmological background. 
Recently, numerical relativity has developed so much that the formation of PBHs can be simulated and the threshold 
of PBH formation can be obtained Shibata and Sasaki @ gave the threshold '00,c — IT — 1-8, where 0n r is 

the threshold value of the peak value 0o of the curvature variable ip and this is consistent with 5c — 0.3 — 0.5 p^[l§|. 
Polnarev and Musco m gave 5c — 0.45 — 0.66 for a radiation fluid, depending on the perturbation profiles, where 5c is 
the threshold value of the averaged density perturbation 5 at horizon entry in the comoving slicing in the lowest-order 
long-wavelength expansion. Musco and Miller [T^ gave 5c for different values of P. Nakama et al. [T^ investigated 
the PBH formation threshold for a much wider class of initial curvature perturbation proHles characterized by five 
parameters. Harada et al. [I^ derived an analytic formula 5c = [3r / (3r -I- 2)] sin^ [ttVP — 1/ (3r — 2)] under a certain 
set of assumptions, which gives 5c — 0.4135 for a radiation fluid. Young et al. |1§| su ggested that there is a signihcant 
environmental effect on 0o,c but not on 5c- Very recently, Tada and Yokoyama and Young and Byrnes [2l| have 
found that the PBH production rate is significantly biased if the statistics of the primordial curvature perturbation 
has small local-type non-Gaussianity. 

To discuss the PBH formation or any other primordial structure formation, it is important to give cosmological 
perturbations which can be generated in the early Universe. On the other hand, the initial data set which will result 
in PBH formation must nonlinearly deviate from the Friedmann-Robertson-Walker (FRW) universe even much before 
the horizon entry of the scale of the perturbations. Such cosmological primordial perturbations are formulated by Lyth 
et al. [ 2 ^ . There have been two approaches so far in which the appropriate initial data of cosmological perturbations 
are prepared and the Einstein equations are numerically solved in spherical symmetry to follow the formation of 
PBHs. The one is initiated by Shibata and Sasaki and the other is by Polnarev and Musco [HI- Although these 
two approaches study the same problem, the direct comparison of their results has not been done and for this reason 
the complete picture of PBH formation has yet been unclear. 

In this paper, without assuming symmetry, we construct cosmological long-wavelength solutions which can deviate 
from the FRW spacetimes with an arbitrarily large amplitude and are naturally generated in inflationary cosmology. 
These solutions can be applied as the initial data for PBH formation and any other primordial structure formation 
in spherical symmetry or in nonspherical situations. Nonspherical cosmological perturbations should be particularly 
important for PBH formation in the stage when the effective equation of state is very soft [23l - l^ . Such stages have 
recently been discussed in some scenarios of inflationary cosmology [23 - [2^ . 

We explicitly show that the solutions in spherical symmetry reduce to the ones in Shibata and Sasaki Q in some 
gauge and to the ones in Polnarev and Musco in another gauge. Thus, we derive the correspondence relation 
between these two approaches. Based on this correspondence relation, we compare the numerical results obtained in 
these two approaches and find that the results are consistent. We analyze three different measures of the amplitude 
of perturbations, 5, 0o a-nd the maximum value of the compaction function Cmax- The threshold values of these three 
quantities depend on the initial density (or curvature) profiles. We analyze how the threshold values depend on the 
profiles and physically understand the dependence using simplified models. We also discuss an environmental or bias 
effect on 0o,c when perturbations are superimposed on perturbations of longer wavelength. 

This paper is organized as follows. In Sec. H, we present the 3-I-I formulation of the Einstein equations and the 
cosmological conformal decomposition. In Sec. HI, we review the long-wavelength scheme. In Sec. IV, we present 
the long-wavelength solutions in different gauge conditions. In Sec. V, we present the gauge conditions adopted in 
Shibata and Sasaki and Polnarev and Musco and find the relation among quasilocal quantities in spherical 
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symmetry. In Sec. VI, we show the equivalence between the long-wavelength solutions obtained in the two approaches 
and derive the correspondence relation. In Sec. VII, we compare and interpret the numerical results obtained in the 
two approaches. In Sec. VIII, we discuss the environmental effect on tpo^c with several simplified models. We use the 
units in which G = c = I and the signature convention (—h -I--I-). 


II. BASIC EQUATIONS 


A. 3-|-l formalism 


We here present the 3-I-I formalism of the Einstein equations according to Nakamura et al. 
in four-dimensional spacetimes is written in the following form: 


|. The line element 


ds^ = —a^dt^ -I- (dx* -|- P’‘dt){dx^ + jd^dt), 


( 2 . 1 ) 


where a, /3® and 7 ^ are the lapse function, shift vector and spatial metric, respectively. The Latin uppercase indices 
run over 1 to 3 and we drop and raise them by jij and its inverse 7 ®-^, respectively, unless otherwise specified. The 
spacetime metric g^i, and its inverse are given by 


= 


A \ 

Pj ) 


( 1 


andff^*^ = 






( 2 . 2 ) 


respectively, where the Greek indices run over 0 to 3. Because of the construction of the inverse matrix, we find 
g = —a^ 7 , where g = det{g^^) and 7 = det( 7 ij). The covariant and contravariant components of the normal unit 
vector to the t =const hypersurface E are given by 


= (-a, 0,0,0) 


and 



(2.3) 


respectively. The projection tensor to E is defined as := g^i, + n^n^. 

The stress-energy tensor for the matter fields is decomposed into the following form: 

= En^^riv + Jp.nv + + Saph^hl, (2.4) 

where E := T^^nf^n'', Ja ■= and Sajs := 

The Einstein equation G^n, = can be written in the following set of equations. The Hamiltonian constraint 

G^'^n^riy = SirT^’^rif^ni, and momentum constraint = 8TTT^^'^n^hi,i reduce to 

n + K^ - = 16 ttE (2.5) 

and 

- V,K = (2.6) 


respectively, where Di and TZ denote the covariant derivative and Ricci scalar with respect to ■jij , respectively, Kij is 
the extrinsic curvature of E defined by 

iL,. := (2.7) 

2 q: 

and K := Note that the semicolon denotes the covariant derivative with respect to 

The evolution equations G^'^h^ih^j = are given by 


Kij^t = a{Tlij+KKij)-2aKiiKj-8TTa 


+ ^%j{E - Sj) 


-V,V,a+(V,/d^)K^,+(V,^^)Kmy+^^V^K,,, ( 2 . 8 ) 


where TZij is the Ricci tensor with respect to "fij on E with TZ = 'y'-^TZij. Equation (12.71) can be rewritten in the 
following form: 


7ii,t — —2aKij + "DjPi -f "EiPj- 


(2.9) 
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The conservation law = 0 is decomposed into = 0 and = 0. Assuming a perfect fluid, 


= {p + p)Ufj.u„+pgfj,„, ( 2 . 10 ) 

where is the four-velocity of the fluid element which is normalized as u^Uf^ = — 1 , they reduce to 

^ jl jm 

= -[^p{v^ + /3')].z + a^pK - ^ (2.11) 

tj + p 

and 

+ {.y/lJmv’’),l = -a^p^m “ ^/7{P + E)a^m H- , , -h (2.12) 

l{p + h) 

respectively, where 


(2.13) 


is the coordinate three-velocity of the fluid element. It is often useful to define the “baryon” or conserved number 
density n and the conservation law = 0 gives 


{y/jau*n)^t + {\/iau*'nv^)^i = 0 . 


(2.14) 


The energy conservation gives 


dn dp 
n p + p 


(2.15) 


The expressions for the four-velocity of the fluid element are rather complicated. In fact, Eq. (I2.13|) implies (Lyth 
et al. [ 2 ^) 


- (/3fc -I- Ufc)(/3'= -I- v'^) ’ 

i t i 
U = U V , 

Ut = -u^[a^ - +Wfc)], 

= u\vi + j3i), 

where Vi = jijvE We can write down A, Ji and Sij as 

E = {p + p)w^ -p, 

2 

Ji = —{p + p){v^+ Pi), 

a 

and 


Sij = Plij + 


E+p' 


(2.16) 

(2.17) 

(2.18) 
(2.19) 


( 2 . 20 ) 

( 2 . 21 ) 


( 2 . 22 ) 


respectively, where 


w.= au^=[l-a [Pk + Vk)[P^-\-v^)] 


(2.23) 


B. Cosmological conformal decomposition 


Here we briefly revisit the cosmological conformal decomposition according to Shibata and Sasaki Q. See also 
Gourgoulhon [^. We decompose the spatial metric 7 ^- into the following form; 


7y = ■0'‘a^(t)7y. 


(2.24) 
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We choose 7 ^ so that 7 = det( 7 ij) is time independent and equal to rj = det(jiij), where rjij is a time-independent 
metric of the flat three-space. The function a(t) is the scale factor of the reference universe. The extrinsic curvature 
is decomposed as 

= A,, + ^K, (2.25) 

where Aij is traceless by definition. We also define A^j as follows: 

.4*^ = i/,-4a-2iy or Aij = ijj'^a^A.j. (2.26) 


Thus, j^'^Aij = 0 by definition. We raise and drop the Latin lowercase indices i,j,k, ... of tilded quantities by 7 ^ 
and 7 ij, respectively. We define Vi and Vi as the covariant derivatives with respect to rjij and 7 ^-, respectively. We 
denote the Laplacians with respect to rjij and 7 ^- as A := rf^ViVj and A := Y^ViVj, respectively. 


Using this decomposition, we can rewrite TZij as follows: 

=n,j+nf^, ( 2 . 27 ) 

where 

2 2 0 2 

Tifj := -—ViVjtfj - —7yA'0 -I- ■^Vi'ipVj'ip - —^ijVk'ipV'^'ip, (2.28) 

:= i[-A 7 ,, + V,V^%, + V.V'^ik, + 2Vk{f^CHj) - 2& (2.29) 

(2.30) 

:= + Vj^a - Vi%). (2.31) 

To derive Eq. (12.291) . we have used the relation 


1 

? 


7"^Vk% = -Vkj = -Vkij = Tt^Vkr]ij = 0. 


The following relations will be useful: 




TV^ - -O' TZ^ = - — 


ViVjip - -^ijAi/j 


6 

■02 


VitfjVjtfj - -^ijV^tpVkii 


The Hamiltonian and momentum constraints are given by 

2 

, - - Zl _L 


nl - AijA^^ + = iGTrE, 


and 


VjAl - -ViK = SttJj, 


3 

respectively. These can be transformed to the following form: 

A0 = - 27r05 

8 8 

2 


AiiA^^ - -K'^ 


V3{rA,,) - -rV,K = SnJi^Z 


The evolution equations become 


{dt - Cp)A,, = ^ [a - (D,D,a - ^DuD^a ] 


1^3 


a^tp'^ 

+a{KA,, - 2A,kA]) - ‘^{Vkp^)A,, - ^ 
{dt - Cp)'4) = + ^{-aK + Vk/3’"), 

{dt - Cp)K = a (AjA^^ + - 


{S.j - , 


(2.32) 

(2.33) 

(2.34) 

(2.35) 

(2.36) 

(2.37) 

(2.38) 


(2.39) 

(2.40) 


DkD^a + 4Tra{E + S^) 


(2.41) 
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where is the Lie derivative along /3® and acts on a scalar field / and a tensor field as follows: 

-C/s/ij = + I3';jkj + P’^jfk^■ 

The definition of the extrinsic curvature yields 


(2.42) 

(2.43) 


{dt - = -2aAij - -%Vkl3'". 

The hydrodynamical equations (12.111) . (I2.12|) and (12.141) are, respectively, written in the form 

{{p + P)w'^ - p}] t + ^ {{P + P)w^ - p] v'-] ; 

= [07V'®a^p('y* + /3*)] ; + ailj^a^pK - + p)(u* + /?*) 

+a-V°aW(p + p)(u'+/3')(u™ + /3’") + , 


(2.44) 


and 


{w-ip^a^ip + p)uj)^t + —^iVvwip^a^ip + p)v'"uj)^k 
\/P 

-atjj^a^pj + wip^a^ip + p) (-waj + Uk^j - 


(wtp^a^n) t H- i^wilPa^nv^) fe = 0. 

’ Vv 


(2.45) 


(2.46) 


(2.47) 


We can define the background Hubble parameter Ht as Hb := a/a. The scale factor a{t) will be that of the flat 
FRW universe in the next section, although the field equations up to here are independent from this choice. 


III. LONG-WAVELENGTH SGHEME 

We here review the long-wavelength scheme based on Lyth et al. [2^ . 


A. Basic assumptions 

To be precise, we focus on some fixed time and put a fictitious parameter e in front of the spatial partial derivative, 
e.g., di —>■ edi. We expand exact solutions in a power series of e, require the field equations at each order and finally 
set e = 1. This scheme is called the gradient expansion. We assume that the spacetime is approximately smooth 
for scales greater than k~^ in terms of the comoving coordinate. This implies that the physical scale of smoothing 
is L = a/k. Since the Hubble length H//^ is the only natural scale of the cosmological evolution, we can make the 
identification e = /L = k/{aHb). Thus, the gradient expansion implies that we assume that the smoothing length 

L is much larger than the Hubble length, i.e., L ^ H//^. 

Based on the above consideration, we make two key assumptions. First, we assume that ip is identically unity 
somewhere in the universe. This makes a{t) be the scale factor for that part of the universe. Second, we assume 
that in the limit e —^ 0, the universe becomes locally homogeneous and isotropic, i.e., an FRW universe, which we 
additionally assume is flat. Thus, the measurable parts of the metric should reduce to the flat FRW one in the 
smoothing scale L which is much larger than the Hubble length This means that there exists a coordinate 

system with which the metric of any local region is written in the following form: 

ds^ = —dt^ + a^pijdx^dx^. (3-1) 

Hence we assume /3® = 0(e), although this is a matter of coordinate choice. As for the spatial metric, a homogeneous 
time-independent ^ij can be transformed away, while a homogeneous time-dependent jij should not exist by the 
present assumption. Since the term of 0{e) in turns out to be decaying, we assume 0{e^). The above 

assumptions are partially justified in the literature in the context of inflationary cosmology [33, 









B. Energy conservation and curvature perturbation 


Let us consider a perfect fluid given by Eq. (12.101) and choose the spatial coordinates so that the worldlines on 
which a;* = const coincide with the worldlines of the fluid elements, which Is called the comoving threading. This 
implies = u^/u* = 0. Then, is given by 


1 


\Jd^ - 


0 , 0 , 0 ) = (-, 0 , 0 , 0 )+ 0 (e") 


= -\Jd^ - 






The expansion of is given by 


• — 7/^ — 


_ 3 

a^ 6 Q, 3 '"' \ y^Q!^ - /3*/3i / OL i>^a 


1 „ / atp^a^ 

Ot 


0 (e 2 


(3.2) 

(3.3) 

(3.4) 


where it should be noted that 7 is time independent. The relation between t and the proper time r along the worldllnes 
of the fluid elements Is given by 


dt 

dr 


1 


The energy conservation law 


Implies 


V'a" - /3'A 

Q = = ^ + {p + p)e 


- + 2 - = + 0 (e 2 ). 


a 'i/i 3 p + p 


(3.5) 


(3.6) 


(3.7) 


If we choose the uniform-density slicing, on which p = const, and assume that the pressure is homogeneous to 0(e), 
'ip/'ip is also homogeneous to 0(e). Since we assume that ip is identically unity at some point, we can conclude that ip 
is time independent to 0 (e), i.e.. 


^ = 0{e^). 


The expansion of is given by 


6n ■■= nP, = 


3 ('0^a),t 


1 


-A 


(3.8) 


(3.9) 


a ip^a aip^a^^ 
and this is related to the trace of the extrinsic curvature through 

= -K. (3.10) 

Since /3* = 0(e), Eqs. (13.41) and (13.91) imply 

0 = 0„ + O(e2). (3.11) 

This shows the equivalence between 9 and to 0(e). For convenience, we introduce the Hubble parameter H as 


3 a 'ip^a 


From Eqs. (13.41) . (13.61) . (13.111) and (I3.12|) . we find 

dp 


^ = _3i7(p + p) + 0(e2). 


(3.12) 


(3.13) 
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C. Equivalence among different time slicings 

We use the following terminology about the slicing conditions. We call the slicing which is orthogonal to the fluid 
worldlines, namely, = u^, the comoving slicing. This is independent of the comoving threading with which the 
worldline of the fluid coincides with that of the constant spatial coordinates. We call the slicing on which the trace 
of the extrinsic curvature is uniform the constant mean curvature (CMC) slicing. Although this is sometimes called 
the uniform-Hubble slicing, we here avoid this terminology because of the ambiguity in the definition of the Hubble 
parameter in inhomogeneous cosmology. We call the slicing with a = 1 the geodesic slicing. 

The comoving slicing implies Ui = 0 and, hence, Ji = 0. This is possible only if is vorticity free because of the 
Frobenius theorem. This condition is physically reasonable in the early universe because the vorticity is not generated 
from vacuum fluctuation during inflation and because the vorticity is conserved for a perfect fluid in any spacetime. 
Therefore, in the comoving slicing, the momentum constraint (j2.36|) implies 

d,H = 0(e3) (3.14) 

and the Hamiltonian constraint (12.351) gives 

^' = YP + 0(e2), (3.15) 

implying 

d,p = 0{e^). (3.16) 

That is 5p = O(e^) and 5H = O(e^) in the comoving slicing, where 5p and 5H are the inhomogeneous parts of p 
and H, respectively. This means that the uniform-density, CMC and comoving slicings coincide to 0(e). This also 
guarantees that tp is time independent to 0(e) in each of those slicings and the time-dependent part appears only 
from O(e^). 

From Eq. (13.131) . we find 


Lp = -3H{p+p) + 0{e'^). (3.17) 

a 

Since p is homogeneous to 0(e), p is also homogeneous to 0{e). Hence, the lapse function is given by 

a = + 0{e^). (3.18) 

P + P 

where A{t) is a function of t. If p is homogeneous to 0{e), which is the case for the barotropic equation of state, a 
is also homogeneous to 0(e) and we can choose a = 1 by rescaling the time coordinate. Thus, the geodesic slicing is 
also equivalent to the uniform-density slicing to 0{e). 


IV. LONG-WAVELENGTH SOLUTIONS 
A. Gradient expansion 

Here we derive the expansion of the physical quantities based on the assumption of /3® = 0(e) and 7 ^ = O(e^). In 
the limit e —?► 0, the metric functions a and ip are assumed to be locally homogeneous. We can use the scaling of time 
so that a = 1 locally in this limit and this suits the metric form given by Eq. dSU). 

From Eq. (12.441) . we find Ay = O(e^). We are still allowed to take general time slicings on which the density is 
uniform to 0{e). As inhomogeneity appears in the mean curvature K = —from O(e^), Ji = O(e^) follows from 
Eq. (I2.38|) . Combining this with Eq. <\2.21L we find Vi + pii = O(e^). Therefore, since ,5* = 0(e), we find u® = 0(e). 
Since 7 ^ = O(e^), we can assign hy = O(e^) using an appropriate coordinate transformation, where hy := 7 y — r/ij. 
In summary, if we assume 

/3® = 0(e), % = 0(e2), (4.1) 

then we deduce 

ip = ip{x^) = 0(6°), u® = 0(e), 5 = O(e^), iy = O(e^), hy = 0{e^), x = 0{e^), k = O(e^), (4.2) 
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c P- Pb , ~ K - Kb 

0 ■■= -, hij := 'yij - riij, x ■= a - 1, k := 


Pb 


n — Hb 


For later convenience, we define 

5n '■= 

Ub l + R 

where we have used Eq. (12.151) and defined 


^ X I f _ P Pb 

0 + 0[e ), Op 


Pb 


Kb 


= fs + 0 ( 0 ), 


R:=« c^= (f 

Pb \dp 


Equations (12.111) . (I2.37L (12.401) . (12.411) and (12.471) in 0(e°) give 


Kb = -3-, 
a 

a\^ Stt 

= -irPb^ 
a J 3 


- = -'^{Pb + iPb), 
a 3 

Pb = -3{pb+pb)-, 
a 

{a^nb),t = 0 . 


This is a complete set of the Eriedmann equations with spatially flat geometry. 

To discuss the next order term in tp, we decompose it into the following form: 

where dt = 0(e°) and ^ = O(e^). It is also useful to note ui = 1 + 0(e®) from Eq. (|2.23l) . 


(4.3) 

(4.4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 

(4.9) 

(4.10) 

(4.11) 


B. Field equations in general gauge 


As we have seen in the previous section, the comoving, CMC, uniform-density and geodesic slicings coincide to 0(e). 
This ensures that we can consistently apply the gradient expansion in any of these four slicings. Here we derive the 
field equations in general gauge which is consistent with the long-wavelength scheme by generalizing the formulation 
by Shibata and Sasaki Q. 

First, Eq. (12.401) results in 4' = 0. This is the starting point of the formulation. Then, Eqs. (12.451) and (I2.47P yield 

(5 -k 6e + Vkv^ - 3HbR{6 -Sp-x-K) = (4.12) 

and 


respectively, where 


+ + = Oie^), 

1 


Vfcu" := 




Vj/ 6 jyl /2 

Equations dHH), (I^TiHl) . (1^371) . (ITTiH) . dm, (I71T1I1) . (I^:M1) yield 
9t[a^(l -I- R)pbUj] = -a^pb [RdjSp -I- (1 -I- R)djx] + 0{e^)., 
6^-3H6(x + k)-V fc/3'= = 0(6^), 

= -27r4'^aVt,((5 - 2k) + O(e^), 

= i(3i? - 1)« - ^(1 + R)x-^{d + 3R6p) + 0{e^), 

dthij = -2Ay + pik'Dj^’" + Vjk'DiP^ - ^rjijVkld^ + 0(e'‘), 


dt-^ij ~t“ ' ^ij — 2iTr4 

,6 Ai 




+ Oie% 


2Hb-i'^VjK = 87r4''>(l + R)pbUj -k 0{e^), 


(4.13) 

(4.14) 

(4.15) 

(4.16) 

(4.17) 

(4.18) 

(4.19) 

(4.20) 

(4.21) 
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respectively, where we have used 

E + S^,=p + 3p + 0{e^), (4.22) 

S., - = Oie^) (4.23) 

following from Eqs. (I2.20I) ~ (I2.22I) . Equations (14.121) . (14.131) and (14.161) give 

S-Sn- SHbRiS -Sp-x-n) = 0{e*), (4.24) 

S - 3HbR{5 - Sp) + 3(1 + R)Hb{x + k) = C>(e^), (4.25) 

(5 - (1 + R)L - 3HbR{6 - Sp) = 0{e‘^), (4.26) 

where it should be noted that V+ v^) = 0{e^). 


C. Long-wavelength solutions 

1. Aij and hij 


We can see that Aij and hij 


are not sensitive to the choice of time slicings to 0{e^). Equation (14.201) admits the 


following solution for A, 




Aij — Pij ^3 


da 


Hb{d) 


0{e% 


where we have put 


Pijix^) := 


\I/4 


4- 


-- I V.Vi'i' - -VtjA'if 


6 




(4.27) 


(4.28) 


and we have omitted the integration constant because it generally gives a decaying mode. 

The expression of hij depends on the choice of the spatial coordinates. For the normal coordinates, where /I® = 0, 
by integrating Eq. (I4.19|) . we find 


hij 


r dd 


r da 

Jo Hb{d) 


+ Cij + O(e^). 


(4.29) 


We drop the constant tensor Cij so that hij —>■ 0 in the limit t —>■ 0 by choosing appropriate spatial coordinates. Of 
course, other coordinate conditions may be useful. For example, in spherical symmetry, we can even have hij = 0 
identically, for which /3® is accordingly determined. This is called the conformally flat spatial coordinates. 

If we assume T = const, we have the background solution from Eq. (1131) as 


Qj = (Xq 



(4.30) 


and we can explicitly write Aij and hij in the following form: 

^ “(3r + 2)(3r-2)^*^' ( 374 ) 

where the bottom one is obtained in the normal coordinates. Aij and hij depend on time as 

Aij 0 ^t^~^ and hij(xt^~^, 


(4.31) 

(4.32) 


(4.33) 


respectively. 


























12 


Z. CMC slicing 


Hereafter we assume F = const for simplicity. Then, Eqs. dH and ( 1 ^ . respectively, yield 

Sn = ^5 + 0{e^), 

5p = S + 0{e'^). 

We present the solutions in the CMC slicing, where k = 0. Equation (|4.17|) is solved for S to give 


S = - 




27r4'®a^P6 


O(e^). 


Defining 


/(x'=) := - 


4 


3 4-5 

and using Eqs. gj]) and (I4.18E we can express S and x as follows: 

2 


= f f—) 

KaHtJ 


S = f 


3r - 2 / 1 


2 


Then, from Eq. (I4.15L we find 


X = 


1 

Uj = —Oi j — 

^ 3T ^ a 


3 r " \aHb 


da 


C 


L^O 


0(e4). 


1 




(4.34) 

(4.35) 

(4.36) 

(4.37) 

(4.38) 

(4.39) 

(4.40) 


Hb{d) \ \aHb^ 

where C is a constant of integration. If C is nonvanishing, there appears a decaying mode for 1 < F < 2. So we 
assume (7 = 0. Thus, we have 

2 


= + /3j) = 


dd 

Hbid) 


(i)' 


O(e^). 


(4.41) 


We can explicitly prove that the above solves the momentum constraint equation (14.211) to 0{e^) using the identity 

P^(4'%) = 4-^^/. (4.42) 

If we choose the normal coordinates, where /3® = 0, Vj is obtained by Eq. (14.4111 as 


1 A 

Vi = - Oi J — 

^ 3F ^ a 


dd 


( 1 


/o Hb{d)\ \aHb 

In this case, ^ is obtained by integrating Eq. (j4.16l) as 


+ 0(e5). 




+ 0(<!-‘), 


where the integration constant is absorbed into 4^. 
Using the background solution (I4.30p . we find 

2 


<5 = / 

K = 0, 

X = - 


1 




aHb, 

3F-2 / 1 y 
3F ^\aHb) 


0{e% 


Uj = u\vj + 13j) = 




3F(3F + 2) \aHb 

3 


1 


O(e^), 


djfa 


(- 


3F(3F + 2) \aHb 




(4.43) 

(4.44) 

(4.45) 

(4.46) 

(4.47) 

(4.48) 

(4.49) 




+ 0(e"), 


(4.50) 
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where the last two equations are obtained in the normal coordinates, 
summarized as follows: 

6 (X , X cx Uj = u*{vj + (3j) oc , 

where the last two equations are obtained in the normal coordinates. 


The time dependence of the perturbations is 
^ Vj , (4.51) 


3. Uniform-density slicing 


The uniform density slicing implies 5 = 0. By the assumption F =const, we find 6n = 6p = 0. From Eqs. (I4.17p . 
(14.251) and (I4.15p . we find 


K 


X 


Uj 



0{e% 



+ 0{e% 


u\vj +/3j) 


1 r r da 

2 ^■'a [Jo Hbid) 



0{e% 


(4.52) 

(4.53) 

(4.54) 


where / is defined by Eq. (14.3711 . We can explicitly show that the momentum constraint (j4.21l) is satisfied to this 
order. Using the background solution (14.301) . we have 


6 

K 


X 




0 , 



O(e^), 





3r + 2^^ -^“ (aiFfc) 


O(e^), 


(4.55) 

(4.56) 

(4.57) 

(4.58) 

(4.59) 


where the last two expressions are obtained in the normal coordinates, and in the bottom equation we have used that 
^ can be shown to be time independent and is absorbed into 4' to O(e^). Therefore, in the uniform-density slicing 
with the normal coordinates, if is time independent to O(e^). 


4- Comoving slicing 


In the comoving slicing, and hence 






Ui = 0. 

(4.60) 

Then, 

Eq. 

(14.151) imnlies 








(r-i)5 + rx = o(e4). 

(4.61) 

Using 

this 

and Eqs. (14.181) and (14.241). 

we find 






dsK -- 

= i(3r-4)«-i5 + 0(e4), 

(4.62) 




dsS-- 

= -SFK-f 3(r- l)5-bO(e'^), 

(4.63) 

where 

we put s = Ina. The matrix 








/(3r-4)/2 -1/2 \ 

1 -3r 3(r-i) J 

(4.64) 
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has eigenvalues (SF — 2) and 3(r — 2)/2 with associated eigenvectors 



respectively. Therefore, the general solution for k and S is given by 


(4.65) 


(4.66) 


where Ci and C 2 are integration constants. For 2/3 < F <2, the first term grows, while the second term decays in 
time. Hence we drop the second term by choosing C 2 = 0. Thus, we have 


« = -^<5 + 0(e"). 

Substituting this into Eq. (I4.17p . we have 

3T + 2-' \aHbJ ’ 


(4.67) 


(4.68) 


where / is defined by Eq. (I4.37p . It should be noted that the time dependence of 5 is consistent with the first term 
on the right-hand side of Eq. (14.661) . Now it is straightforward to see that k, x and ^ are given by 


K = 


X = 


- f 

sr + 2-' 
3(r-i) 
3r + 2 




0(e% 


Uj = 0 , 

Vj = 0 , 

C =- ^ -/ 

^ 2(3F-H2)-' 



(4.69) 

(4.70) 

(4.71) 

(4.72) 

(4.73) 


where the last two expressions are obtained in the normal coordinates. We can explicitly show that the momentum 
constraint ()4.21l) is satisfied to this order. The combination of the comoving slicing and comoving threading, i.e., 
Ui = Pi = 0, is called the comoving gauge. 


5. Geodesic slicing 

In the geodesic slicing, it is straightforward to obtain 




9r-4'' \aHb 


1 


0{e% 


3r-2^ 
K = - rf 


9F-4-' \aHb 


1 




X = 0, 

6 (r-i) 

Uj = — ——Ojja 


(_L\ 

i9T-4){3T + 2y^'’^\aHbJ 


6(r-i) , 

{9T-'l){3T + 2) ^^^\aHb 
2 


Vn = - 


e = - 


1 


/ 


2(9F-4)'' \aHb 


1 


+ 0(e4), 


0{e% 


O(e^), 


(4.74) 

(4.75) 

(4.76) 

(4.77) 

(4.78) 

(4.79) 


where the last two expressions are obtained in the normal coordinates. The gauge condition of a = 1 and /3® = 0, i.e., 
the geodesic slicing with the normal coordinates, is called the synchronous gauge. 
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6. Relations among the different time slicings 

In Sec. imi we have seen that the slicings in the long-wavelength scheme coincide to 0{e). This means that the 
zeroth-order curvature variable d' = should be common, i.e., 

^'CMC = 'I'UD = ^*0 = ^*0 (4.80) 

for the same spacetime, where CMC, UD, C and G stand for the CMC slicing, uniform-density slicing, comoving 
slicing, and geodesic slicing, respectively. Using the results obtained in Sec. II V Cl we can find the relations among 
the perturbations in the different slicings. Irrespective of the choice of /3*, we can conclude 


= Atj^c + 0{e^) = A^j,G + 


(4.81) 


and 


5ud = 0, = ttf;—7<5cmc, 


3r-b2 


9r-4 


2 2(3r — 2) 

«:CMC = 0,KC = 2 + 0(6"^), kq = ^ «ud 


Xc = 


9r(r — 1) 3r 

XCMC + 0(€^), xuD = — ^vi^CMC + 0(e^),XG = 0, 


2(3r-2)^ 


(3r-2)(3r-b2) 

n 3r 4 , 9r(r-i) 

Uj C = 0, Uj UD = CMC + C(e },Uj G = - QP _ 4 cmc- 

In the normal coordinates, where /3* = 0, we can also conclude the following relations: 

hij,CMC = hij^cD + O(e^) = hij^G + O(e^) = + 0(6"^), 

3r 4 9r(r — i) 

Vj UD =-CMC + 0{€^),Vj c = 0,Vj G = - gp _ ^ Vj CMC, 

3r 3r 

= OT^ , o ^CMC +0(e^),^uD = 0(e^),^G = tg;-t^Cmc- 


3r-b2 


9r-4 


(4.82) 

(4.83) 

(4.84) 

(4.85) 

(4.86) 

(4.87) 

(4.88) 


V. SPHERICALLY SYMMETRIC SPACETIMES 

One of the important motivations of the following sections is to compare the results of numerical simulations in 
two very different approaches to PBH formation in spherical symmetry. The one uses the CMC slicing with the 
conformally flat spatial coordinates, while the other uses the comoving slicing with the comoving threading, which 
is called the Misner-Sharp formulation of the Einstein equation in spherical symmetry with a perfect fluid. For this 
reason, we review the two formulations of the Einstein equations in spherical symmetry. 


A. Gauge conditions in the two approaches 

Shibata and Sasaki Q adopt the conformally flat spatial coordinates. The line element is given by [sj 

ds^ = —{a^ — f3‘^r'^)dt^ + 2ifj^a^Prdrdt + {dr'^ + r^dVi^), (5.1) 

where = dO^ + sin^ Odcj)^ is the line element on the unit two-sphere. Thence, we have /?’' = r/3 and Xij = Vij- 
The spatial metric is conformally flat and this is also a minimal distortion gauge because jij = 0. As for the slicing 
condition, Shibata and Sasaki Q adopt the CMC slicing on which K = —3Hb{t) and formulated the initial value 
problem of the Einstein equation. The combination of the conformally flat spatial coordinates and the CMC slicing 
fixes the gauge. Since the full set of field equations can be derived from the equations in Sec. Ill Bl and are given in 
Shibata and Sasaki Q, we do not repeat them here. 

If we adopt the comoving slicing with the comoving threading, we have a compact set of field equations and this is 
called the Misner-Sharp formulation [s^. The line element is given by 


ds'^ = —a^dt^ -I- b'^dr'^ + R^dVt^. 


(5.2) 
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The following equations derive from Eqs. (j2.5l) . (12.6L (I2.7L (12.81) . (I2.11L (j2.12|) . (I2.20|) and (j2.21l) with 

M' = AirpR^R', 

M = —AnpR^R, 


R' = R— + R', 
a 0 

V = -(/0 + p) — 

Q 




1 - 


62 a2 


where the prime denotes the partial derivative with respect to r and M is called the Misner-Sharp mass 
formulation has been adopted for PBH studies by many authors including Polnarev and Musco |ll| . 


= z;* = 0: 

(5.3) 

(5.4) 

(5.5) 

(5.6) 

(5.7) 
This 


B. Areal radius, mass excess and compaction function 


In spherical symmetry, the areal radius is defined by i? = where A is the area of the 2-sphere with 

constant t and r. We define the Kodama vector by = e^^dsR, where cab is the totally antisymmetric tensor 
in the two-dimensional manifold charted by t and r and the Latin uppercase indices run over 0 and 1, and cab is 
given by eab = —Geab with the Levi-Civita symbol eab and the determinant G of the two-dimensional metric 
Gab- We raise the indices of cab as = cacG^^ and = G^^e(?. We trivially extend to as a vector 
field in the four-dimensional manifold. Since = TjfK'^ is a conserved current, 


M ■=- 


= 




S^ay/^dx^ 


(5.8) 


is a conserved mass, where in the last expression E is chosen to the constant t hypersurface and the integration is 
done within a ball on E. M is called the Kodama mass. This is equivalent to the Misner-Sharp mass in the present 
setting. 

Shibata and Sasaki Q introduce the notion of a compaction function. Here we define it in more general settings. 
To this end we first define an excess 6M in the Kodama mass by the difference between the Kodama masses of the 
two spheres of same t and areal radius R (not r) in two different spherically symmetric spacetimes. This definition 
of SM is covariant with respect to the choice of spatial coordinates but does depend on the slicing. We use the flat 
FRW spacetime as a reference spacetime to define the mass excess. 

In an analogy with asymptotic flatness, if initial data approach those of the flat FRW one so fast that the mass 
excess has a finite limit as we take r oo, we call such data asymptotically flat FRW. For asymptotically flat FRW 
data, if the mass excess approaches zero as we take r —> oo, we shall call such initial data compensated] otherwise, we 
shall call it as uncompensated. It is not clear which model is more physically realistic, compensated or uncompensated. 
Clearly, uncompensated models are more generic than the compensated ones. Physics should be as generic as possible 
within the framework of asymptotically flat FRW data. On the other hand, uncompensated models might be regarded 
as perturbations of infinitely long wavelength. 

We define a compaction function C(t,r) by the ratio of the mass excess to the common areal radius of the two 
spheres, i.e.. 


C{t,r) 


dM{t, r) 
R{t,r) 


(5.9) 


By definition, C(t, r) = 0 for the flat FRW spacetime. For asymptotically flat FRW data, we find C oc I/i? as we fix t 
and take the limit r —>■ oo. In the following we express the compaction function in terms of the density perturbation 
first in the Misner-Sharp formulation and later in more general slicings. 


1. Comoving slicing 


In the Misner-Sharp formulation, for which the line element is given in the form of Fq. dEa), the components 
of the Kodama vector are given by 


A* = - 


r 


R 

ab’ 


ab 


K 


= 0 . 


(5.10) 
















17 


Then, we find 




where we have used = 0 in this coordinate system, and hence 

M(t,r)=A'K f dxpR‘^{t,x)R'{t,x). 
Jo 


(5.11) 


(5.12) 


Note that this expression is fully valid even in the general spherically symmetric spacetime and shows the equivalence 
between the Kodama mass and the Misner-Sharp mass. We denote the Kodama mass in the flat FRW spacetime by 
Mp-, which is given by 


Mp{t,r) = 4:'Kpba^ 



Then the mass excess is defined by 


SM{t, r) := M{t,r) — Mp{t, r), 


(5.13) 


(5.14) 


where f in the FRW spacetime gives the same areal radius R{t, r) of the sphere of r on the constant t hypersurface 
in the perturbed spacetime, namely. 


a{t)f = R{t,r). (5.15) 

Since E is the constant t hypersurface in both unperturbed and perturbed spacetimes, we find 

a{t)df = R'(t,r)dr, (5.16) 


therefore. 


MF{t,f) = Appb f dxR^{t,x)R'{t,x). 

Jo 

Thus, we find the following exact relation: 

SM{t,r) = Airpb f dxR^{t,x)R'{t,x)S{t,x), 
Jo 


(5.17) 


(5.18) 


where 


S.-P_P1, (5.19) 

Pb 

This is independent of the choice of the spatial coordinates. Thus, we have established the theorem that in spherical 
symmetry and in the comoving slicing, the density perturbation integrated over the perturbed spatial geometry exactly 
coincides with the mass excess. We can also have the exact expression for the compaction function as follows: 

C = ^S{HbRf, (5.20) 


where the averaged density perturbation d is defined as 

^ J SdE in dxR^{t,x)R'{t,x)S{t,x) 

’ J dT, 47r Jq dxR’^(t,x)R'{t,x) 


(5.21) 


with the spatial integration in the perturbed metric and the Friedmann equation (HTl) is used. 
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2. General time slicing 


Next we will see the situation in the general time slicing with the conformally flat spatial coordinates, for which 
the line element is given in the form of Eq. (ED. Then, R is given by 


R = 


and the Kodama mass is given by 


from Eq. (1^ . Since are given by 


M{t,r) = 4:TT f dxx^a^aTp^Tj^K^ 
Jo 


K^= -= K<^=K^ = Q^ 

aip"^ aip‘‘a 


M is expressed as 


In the flat ERW spacetime, this reduces to 


(5.22) 


(5.23) 


(5.24) 


pT 

M (t, r) = 47ra^ J dxx^ilj'^ [(p + p)u*ut + p] {ip'^x)' + {p + —(y’^a),t| ■ (5.25) 


MFit,r) = Aira^Pb f dxx^. (5.26) 

Jo 

We have defined the mass excess 6M{t, r) as the difference between the Kodama masses enclosed in two spheres with 
the same areal radius in the perturbed spacetime and the FRW spacetime, i.e., 

6M{t,r) := M{t,r) - Mpit^ilPr). (5.27) 

Assuming the long-wavelength solutions, we find 


M{t,r) = Ana^Pb J dxx^ {1 + 6)tj;^ ^ + 0{€^), 


(5.28) 


while the integral in Mp{t^il}'^r) can be transformed as follows: 

ry^t{r) 


dxx"^ = J dxx^ = J dy^^^ipl = J dyy^tp^ ^1 -|- y (5.29) 


where ifitir) := 'tp'^(t,r)r. Therefore, we find 


5M = Ana^pb J dxx^ip^ ^1 -|- ^ + 0{e^). 


(5.30) 


Thus, we have established the theorem that the density perturbation integrated over the perturbed spacetime coin¬ 
cides with the mass excess to O(e^) in any time slicing which is compatible with the long-wavelength scheme. The 
compaction function satisfies 


C = \AHbRY + 0{e^). 


(5.31) 


We can confirm that in the comoving slicing, where Ur = 0, Eqs. (15.301) and (I5.3ip hold without the term of O(e^) on 
the right-hand side. 
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VI. SPHERICALLY SYMMETRIC COSMOLOGICAL PERTURBATIONS 

A. Long-wavelength solutions 

Now that we have the long-wavelength solutions without symmetry in Sec. IIVI it is straightforward to apply those to 
spherical symmetry. For convenience, we use the spherical coordinates in which the flat metric rjij takes the following 
form: 


rjijdx^dx^ = -f r‘^{d9‘^ + sin^ Odcji^). 


( 6 . 1 ) 


In the CMC slicing, the solution is given by Eqs. (14.311) . (I4.32p . (I4.45I) - (I4.50I) . where the expressions for vj and 


hij hold only for the normal coordinates and 


vp = vl/(r), 

4 1 1 d 


/ = - 


3 dt® dr 


,d^ 

dr 


Pij = 


1 


2 






J 




) 


( 6 . 2 ) 

(6.3) 

(6.4) 


In the comoving slicing, the solution is given by Eqs. (I4.31L (I4.32p . (I4.68l) - (l4.73p . where the expressions for Vj 
and hij are valid in the normal coordinates and / and pij are given by Eqs. (I6.2I) - (I6.4I) . 


B. Asymptotic quasihomogeneous solutions 


Polnarev and Musco m presented asymptotic quasihomogeneous solutions as cosmological nonlinear perturbations 
based on the Misner-Sharp formulation. We briefly review these solutions below. 

We choose the flat FRW spacetime as the reference spacetime, which is given by Eq. (j5.2p with 

a = 1, b = a{t), R = a{t)r = p = Ph{t), M = Mj, := —U = dr = Ht = 

6 a 

where U := R/a. The line element is written in the form 

ds^ = —dr + a^{t)[dr^ + r‘^{d6‘^ + sin^ Odc/T)]. 

We define tq as the comoving scale under consideration. Then, Rq = a{t)rQ is the corresponding physical scale. We 
define 

HbRoJ \aroJ \Hbaro 

Note that e is time dependent and 

- = (3r - 2)Hb (6.6) 

e 

implies that e increases in time if F > 2/3. The correspondence between the Shibata-Sasaki e and Polnarev-Musco e 
is thus given by ~ e, although the perturbation scheme is somewhat different; Polnarev and Musco [ll| introduce 
e as a small time-dependent function, while Shibata and Sasaki Q introduce e as a constant order parameter which 
controls the order of spatial gradients at the fixed time and which is finally taken to be unity after the expansion. 

As for nonlinear cosmological fluctuations, we assume that the metric approaches 




ds^ 


—dr + a? it) 


dr"^ 

1 - K{r)r^ 


+ r^{d0^ +811)^ edcjT) 


in the limit e —> 0 and expand a, b, R, p, U and M as 

a = I -b ea, b = 


— K{r)r- 


yi -|- £&), R — I?b(I + £.R)j 


47r 


p = pbil + ep), U = HbRil + eU),M = y Pbi?^(l + eM), 


(6.7) 
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noting that K{r) = O(e^). We expand the Einstein equations (I5.3I) - (I5.7I) in a power series of e. We can find that the 
first-order functions can be written in terms of -K(r) through Eq. (15.7|) . The concrete expressions are as follows: 


a = 

b = 

R = 

P = 
M = 
U = 


-h 0{e), 


3(r- 1) {r^K{r)y 
3T + 2 ^ 

3(r-l) \{r^K{r)y 

{3r-2){ST + 2f 

1 


(3r-2)(3r-h2) 
3r {r^K{r)y 
3T + 2 ^ 


3r^ 
(E-l) 


0(e), 


{r^K{r)y 


K{r 


+ 0 (e), 


0(e), 


3r 


3T + 2 
1 


3r + 2 


K{r)rl+0{e), 

K{r)rl+0{e), 


( 6 . 8 ) 

(6.9) 

( 6 . 10 ) 

( 6 . 11 ) 

( 6 . 12 ) 

(6.13) 


where we have dropped a decaying mode. 

The density perturbation 6{t,r) averaged over the region inside the sphere of the radius r is calculated as 


= e^^-^K(r)rl + 0{e^). (6.14) 

Expanding S{t,r) as S{t,r) = e6{t^r)^ we find 

- 

= ^^^^ K{r)rl + 0{e), (6.15) 


which is time independent in the limit e —>■ 0. This coincides with the mass perturbation M(t,r) to O(e^), i.e., 
S = M + 0(e). This is consistent with the full-order equation (15.201) noting Mb = H^Ri/2 and Rb = ar. 

For variable equations of state and for more general modes, see Polnarev and Musco [ij. For higher-order solutions 
with respect to e, see Polnarev et al. M- 


C. Equivalence of the two solutions 


As we have seen in Sec. 1111 Cl the CMC slicing and comoving slicing coincide with each other to 0(e). This means 
that the spatial metrics dY? in these slicings are identical to 0(e). In the conformally flat coordinates, the spatial 
metric is given by 

dY^ = [dzj'^ + w^{d9‘^ + sin^ 0d(/)^)] -I- O(e^), 

while in the areal radial coordinates, it is given by 


dY^ 


a^{t) 


dr'^ 

1 - A:(r)r 2 


-f r^{de^ -f sin^ Odcj?) 




where and hereafter we denote the radial coordinate in the conformally flat spatial coordinates with w to distinguish 
between the two coordinates. The former and the latter are adopted by Shibata and Sasaki Q and Polnarev and 
Musco m , respectively. The coincidence of the two metrics implies the following relations: 


'i’{vj)'^dzu = 


dr 


\J\ — K{r)r^ 


'i’{wyw = r. 

The above expression can be inverted for r and K(r) in terms of w and ’^{w) as 

' r = 


K{r)r^ = 1 - 1 -k 2 


d'^{vj) 


'I'(r:7) dw 


(6.16) 


(6.17) 
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Using Eqs. (16.161) and (I6.17|) . we can prove the following relation; 

1 1 d{r^K{r)) Aid 
3 r'^ dr 3 dw 

Let us prove this identity. From Eq. (I6.17|l . we find 


,dAi 

dw 


\Jl — K{r)r^ = 1 + 2vj 


1 dAJ 
'h dw' 


Together with Eq. (I6.16|) . we find 


On the other hand, we hnd 


1 d 


37-2 dr ('h + 2w^) dw' 


r^K{r) = -Aw ^-— I T + w -— ) 
dw \ dw j 


Differentiating the right-hand side with respect to w, we find 


_d_ 

dw 


, dT 


dw 




w -— dt -I- w -— 


dw 


dT 


= 1 -f 2w— ] ^ 

dw / dw V dw 


, dT 


From the above equations, we finally obtain Eq. (I6.18E 
From Eqs. (16.31) and (j6.18|l . we find a useful relation 

/ = 

ar 

Equations (16.161) can be integrated to give w and dt (w) in terms of r and K (r) as 


w = r exp 

= exp 

where for simplicity we have assumed 


dx 


1 


oo ^ \ 

I f'' (^( _ l_ 


- 1 


X y a/1 — K{x)x‘^ 


- 1 


lim K{r)r^ = 0, lim 'I'(-n7) = 1. 


(6.18) 


(6.19) 


( 6 . 20 ) 


( 6 . 21 ) 


( 6 . 22 ) 


In spite of the above equivalence, it should be noted that the turning point where K{r)r^ = 1, which is a coordinate 
singularity in the areal radial coordinates, can be overcome in the conformally flat coordinates as Kopp et al. [s^ 
point out. 

The correspondence between the Polnarev-Musco variables and Shibata-Sasaki variables is given by 

5 = ep + 0{e^), (6.23) 

X = ea + 0{e^), (6.24) 


= -e\[3a - (3r - 2){3R + rR' + b)] + 

O 


(6.25) 


Using the relation ()6.20|) . we can see that the Polnarev-Musco asymptotic quasihomogeneous solutions given by 
Eqs. (I6.8I) - (I6.11[) are equivalent to the long-wavelength solutions in the comoving slicing given by Eqs. (I4.68|) - (I4.70D . 


D. Correspondence relation between the two solutions 

The difference between the Shibata-Sasaki and Polnarev-Musco formulations is in the choice of time slicing and 
spatial coordinates. In Eq. (14.821) . we establish the correspondence relation 


dc = 


3r 


3r-f 2 


<5cmc + 0(e^)- 


(6.26) 
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The prefactor on the right-hand side is due to the difference in time slicing. Note that Eq. (I6.26|) implies 




3r 


ri^CMC + 


3r-H2 

where (5cmc = e<^CMC- In particular, we find the useful relation from Eqs. (I6.15|) and (16.271) . 

^CMC = K{r)rQ + O(e^). 


(6.27) 


(6.28) 


It is also interesting to derive the expressions for physical quantities in terms of both the Shibata-Sasaki variables 
and Polnarev-Musco variables. Shibata and Sasaki measure the amplitude of the perturbation by 'I'(O). This can 
be expressed in terms of the Polnarev-Musco variables as follows: 


^<(0) = exp 


1 

rdr( 

2 

Jo “1 


1 




- 1 


Since 


SM = Ana^pb J -I- ^ + 0{e^), 


in both slicings, Eq. ()6.26l) implies 


SMc = 


3r 


-6Mcmc + O(e^). 


3r-b2 

For the CMC slicing, substituting Eqs. (14.451) and (16.31) into Eq. ()6.30l) . we find 

Aira^p, 


(6.29) 


(6.30) 


(6.31) 


SMcmc = dx{-A)x (24-' -b xAf") (4- -b 2x4'') -b 0{e^) 

= -2007^4'' (4' -b tx74'') -b O(e^). 


(6.32) 


Therefore, the mass excess uniformly increases in proportion to the scale factor. The compaction function C is then 
calculated as follows: 


CcMC = 


i5Mcmc _ 1 

w^'^a 2 


1-1 + 2 


w 


4' dw ) 


+ 0{e^). 


This is time independent to O(e^). From Eq. (16.311) . we have 

3r 


Cc = 


CcMC + 0(c)- 


3r + 2 

From Eqs. (16.17^ and (j6.33|) . we find a very simple relation 

Ccyic = \K{ry 

Hence, from Eqs. (|6.15|) and (j6.27p . the relation between 5 and the compaction function C is given by 

2 


C = i|(r) 


ro 


+ 0{e^) 


and, in particular, 


C{t,ro) = ^5(?'o) + O(e^) 


(6.33) 


(6.34) 


(6.35) 


(6.36) 


(6.37) 


in any slicing. In fact, we can see that Eqs. (16.361) and (16.371) directly follow from Eqs. ()5.20l) and (15.311) . 

Before moving onto the comparison of the numerical results, we show how the asymptotic condition restricts the 
functions 4'(n7) and K(r). If we assume 4'(n7) —>• 1 and K(r)r^ —>■ 0 in the asymptotic region, from Eq. (16.331) . the 
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asymptotic flat FRW condition implies — 1 = O , while the compensation condition implies faster falloff. If 

we assume vj; = l -|- C-^ji^w) + O (l/ro^) , we find 

lim SMqmc = aCjjj, Ccmc ~ (6.38) 

r—)-oo ZU 

The compensated model corresponds to C^p = 0. In terms of K{r), we find the asymptotic flat FRW condition implies 
K{r) = O (l/r^), while the compensation condition implies faster falloff. If we assume K{r) = + O (l/r"^), 

we find 

lim SMcmc = clCk, Ccmc ~ ■ (6.39) 

r—¥oo T 

The compensated model corresponds to Ck = 0. 


VII. NUMERICAL RESULTS 


A. Setup of numerical simulations 


In the decelerated expansion, such as in the radiation-dominated era, the Hubble horizon expands in terms of 
the comoving scale. In this phase, primordial cosmological perturbations which were on superhorizon scale enter 
within the Hubble scale, which is called horizon entry. Much before the horizon entry, such primordial cosmological 
perturbations are naturally described by the long-wavelength solutions. The argument of Jeans instability strongly 
suggests that the density perturbation which collapses to a black hole must be of order unity at horizon entry because 
the Jeans scale is comparable with the Hubble scale if F — 1 = 0(1). Since J(ro) approximately gives the averaged 
density perturbation at horizon entry, the discussion in Sec. IVI PI implies that the compaction function C and the 
curvature variable must have been perturbed by order of unity even much before the horizon entry. Thus, it is 
natural to prepare the nonlinear long-wavelength solutions as initial data sets for PBH formation at the moment much 
before the horizon entry. 

As we have seen, there are two major approaches, the one is adopted by Shibata and Sasaki Q and the other is 
adopted by Polnarev and Musco [ll|. They are different not only in time slicing and spatial coordinates but also in 
amplitude measures and in initial curvature profiles, which lead to the complexity in the comparison of the numerical 
results between these two approaches. Since we have already discussed the formulations, we will focus on the difference 
in initial curvature profiles. Shibata and Sasaki Q determine the initial data by performing iteration to solve the 
Hamiltonian and momentum constraints with the conformally flat spatial coordinates so that the density profile is 
given by 

V’^'^CMC = Cs 

with the boundary condition 

V' = 1 + ^ + (^-2) 

as tu —>• oo, while Jcmc and Uj cmc are those for the long-wavelength solutions, i.e., Eqs. (14.451) and (14.481) . wq is 
fixed and then a and Cf, parametrize the prohles. Note that the above model is approximately compensated because 
the overdensity in the first term and the underdensity in the second term in the square brackets on the right-hand 
side of Eq. (EH cancel out if they are integrated over the whole flat space with ■0 = 1. However, since the space is not 
flat or 0 is not identically unity, this model is generally uncompensated. Another complexity comes from the different 
amplitude measures. Shibata and Sasaki Q adopt 41(0) and Ccmc, max as amplitude measures, where Ccmc, max is the 
maximum value of Ccmc. 

Polnarev and Musco [1^ and Musco and Miller [l^ give Gaussian-type profiles and top-hat-type profiles. In both 
cases, the profiles fall off much faster than as r —>■ oo and, hence, their prohles correspond to exactly compensated 
models. In the Gaussian type, A"(r) is given by 


VJ 

exp [ 

0 


'LkJ r 


exp 


G^VJr 


t 


(7.1) 


K{r) 



(7.3) 


which is parametrized by A and a. They adopt Jc.o •= ^c(^o) as an amplitude measure, where tq is specihed as the 
smallest positive root of p{r) = 0. We will abbreviate 4'(0), Ccmc, max and Jc,o as 0o, Cmax and 5, respectively. 
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Hereafter we only keep the lowest order terms and neglect higher order terms in terms of e of the long-wavelength 
solutions as initial data sets because this truncation is numerically justified if e is sufficiently small in Polnarev and 
Musco [ll|. Equation (16.1111 implies that tq is a positive root of the equation {r^K{r))' = 0, while Eq. (16.3511 implies 
that ri is a positive root of the equation {r‘^K(r))' = 0, where r = ri is the radius at which Ccmc takes a maximum. 
If we assume a top-hat shape for K{r), we can see tq = ri and, hence, 

QF -I- 2 - 

C„.ax = -^5. (7.4) 

61 

In Sec. IVIII C[ we will see the top-hat curvature model in more detail. For an exact Gaussian model, where 

7!:(r) = exp , (7.5) 

we find rg = and ri = y/2A and hence 

Cmax = ^ VeS. (7.6) 

Note that the above relation is justified only for the exact Gaussian profile for K{r). There is no one-to-one corre¬ 
spondence between 'ipQ, Cmax and S unless we specify the profile. 


B. Comparison of the numerical results 


To test the consistency between the results of Shibata and Sasaki Q and the results of Polnarev and Musco m 
and Musco and Miller [^, we will reproduce the former simulation with the latter formulation. The matter field is 
assumed to be a radiation fluid; i.e., P = 4/3 in this subsection. 

To make discussion clearer, we first integrate 




which is obtained by Eqs. (|4.45p and (|6.2p . with the source term 

^2 ' 


= Cs 


w 

exp ( - — 
0 


tUn 


.-3 . 


VO 


-G exp 

\VOo 


(7.7) 


(7.8) 


and the boundary condition 


^ = 1 + ^ + C>(tz7-2) 


(7.9) 


as tn —>■ oo. Then, we determine K{r) through Eqs. (I6.17|) . Equivalently, one may define h{r) := zu{r)/r and solve 


2h 


1-f 


rh' 


2h' iCs 


h{rh) 


'3 


exp(—— cr exp — 


2 ;, 2 \ n 


r^h 


(7.10) 


with the boundary condition h = \ — C-ifj /r + 0(r ^)asr—>-oo. It would be useful to note that h (0) = 4^ (0) ^. Then, 
K{r) can be calculated by 


K{r) 


l( h{rf \ 
r-2 (r/i(r))'2 ) ' 


(7.11) 


The obtained K[r) generates asymptotic quasihomogeneous solutions given in Sec. I VI Bl and we implement numerical 
simulations based on the Misner-Sharp formulation adopting those solutions as initial data sets. Several examples of 
the initial data sets for both ^'(ccj) and K(r) at the black hole threshold are plotted in Fig. [T] for different values of 
cr. The initial data sets constructed here are identical with those in Shibata and Sasaki to 0{e). 

The details of the numerical code are described in Nakama et al. [l^. The threshold models for black hole formation 
for different values of a are summarized in Table U and Fig. [21 where the scales are chosen so that a = aft^, o/ = 1 
and 070 = 1 ■ We can see that the obtained numerical results are fairly consistent with Shibata and Sasaki ’s results 
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K(r) 



FIG. 1. Initial profiles at black hole threshold for the model given by Eq. (iTlSl) with the different values of a in terms of (ro) 
(left panel) and K{r) (right panel). 


TABLE I. Black hole formation thresholds for the model given by Eq. (iTlSl) parametrized by a and Cs- 


a 

1.25 

1.5 

2 

3 

5 

8 

12 

20 

Cs 

31.0 

18.9 

13.2 

10.9 

10.4 

10.9 

11.4 

12.1 

Cv, 

-0.22 

-0.24 

-0.28 

-0.33 

-0.41 

-0.50 

-0.57 

-0.63 

00 

1.40 

1.41 

1.42 

1.46 

1.52 

1.59 

1.63 

1.69 

Cmax 

0.43 

0.42 

0.42 

0.41 

0.39 

0.39 

0.38 

0.38 

6 

0.56 

0.55 

0.54 

0.52 

0.49 

0.46 

0.44 

0.42 


by comparing their Figs. 2 and 7 with our Table HI Note that Shibata and Sasaki Q estimate that the threshold 
value of '00 is — 1.79 in the limit cr —oo. Here we have also calculated and listed the values for and S. We can 
see that is negative for these threshold models, implying that the models are overcompensated. In fact, we have 
numerically found that is negative also for all nonthreshold models we have calculated for 1.25 < a < 20. The 
threshold values of 6 are comparable with those obtained by Polnarev and Musco m and Musco and Miller [l^ . 

It is also interesting to calculate 0o and Cmax for the results of Polnarev and Musco m and Musco and Miller [I^ . 



a 

FIG. 2. The threshold values of tpo, Cmax and 5 for the model given by Eq. (17.81) for different values of a. 
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In fact, this is possible without reproducing the numerical simulations because the initial profiles in Polnarev-Musco 
and Musco-Miller simulations are given explicitly in terms of elementary functions. The results for the Gaussian-type 
proHles are summarized in Table in and Fig. El 


TABLE II. Black hole formation thresholds for the model given by Eq. (El parametrized by a and A 


a 

0.00 

0.25 

0.50 

0.75 

1.00 

2.00 

5.00 

A 

1.01 

0.89 

0.80 

0.73 

0.68 

0.54 

0.36 

ro 

1.75 

1.70 

1.61 

1.51 

1.42 

1.15 

0.80 

ri 

1.43 

1.40 

1.35 

1.29 

1.22 

1.01 

0.71 

Ipo 

1.61 

1.59 

1.57 

1.57 

1.55 

1.53 

1.50 

^max 

0.37 

0.37 

0.38 

0.38 

0.39 

0.39 

0.40 

s 

0.45 

0.46 

0.47 

0.47 

0.48 

0.49 

0.50 



a 


FIG. 3. The threshold values of ipo, C max 


and 5 for the model given by Eq. (El with different values of a. 


C. Threshold values of <5, Cmax and ipo 

Hereafter, we denote the threshold values of S, Cmax and ipo as 6c, Cmax,c and ipo^c, respectively. In Table |T] and 
Fig. n we can see that the larger the tr is, the smaller the Sc becomes. This tendency is shared by Cmax.c- This 
can be interpreted as follows, tuo gives the scale of the overdense region, while atuo gives that of the underdense 
region surrounding the overdense region and gives the ratio of the density perturbation of scale atug to the that 

of scale zuo- If cr is not so larger than unity, the overdensity must be approximately compensated by the relatively 
narrow underdense layer, resulting in a sharp transition to the background FRW universe (see Fig. Ej). This enhances 
the effect of the pressure gradient and suppresses the gravitational collapse, resulting in the larger threshold values 
Sc and Cmax,c- If O' is much larger than unity, the underdense region spreads to large distances and hence the density 
there is only slightly lower than the FRW spacetime. Such a configuration minimizes the effect of pressure gradient 
force. This is the reason why Sc and Cmax,c take minimum values ~ 0.42 and ~ 0.38, respectively, for a = 20. 

The same tendency can be seen also in the Gaussian-type profiles for K(r) as seen in Table IH] and Fig. El The 
larger the a is, the larger the Sc becomes. Figure 7 in Musco and Miller [l^ shows that the larger the a is, the larger 
the density gradient and hence the pressure gradient become. The minimum value 0.45 of Sc is realized for the pure 
Gaussian profile of K(r), i.e., a = 0. The factor ^/e/2 ~ 0.82 in Eq. (17.61) explains the ratio of Cmax,c — 0.37 to 
Sc — 0.45. As for the sharpest transition model, which is an approximately top-hat curvature model, Polnarev and 
Musco [ll| even get Sc — 0.66, which is nearly the possible maximum value 2/3 of 5 for the density perturbation 
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On the other hand, the analytic threshold formula obtained in Harada et al. [T^ yields Sc — 0.4135 for a radiation 
fluid. One of the key assumptions to derive this analytic formula is that the effects of pressure gradient force due to 
the transition from the overdensity to the underdense layer can be neglected. From the above consideration, we can 
say that the analytic formula approximately gives the lowest value of Sc, which is realized if the transition between 
the overdense region and the FRW universe is sufficiently smooth. 

For the equation of state p = (F — l)p, the above discussion about the threshold of PBH formation in spherical 
symmetry is summarized as follows: 


^c,min ^ Sc ^c.max 


in the comoving slicing, where 


5 


c,min 


max 


3r . 2 /ttVf^n 
3F + 2 V 3r-2 ) ’ 
3r 

3F + 2' 


(7.12) 


(7.13) 

(7.14) 


The smoother the transition from the overdensity to the FRW universe is, the smaller the Sc becomes. The minimum 
value is realized if the transition is sufficiently smooth, while the maximum is realized if the transition is sufficiently 
sharp. 

As seen in Tables U and im it does not seem to affect Sc and Cmax.c so much whether the model is compensated or 
not. This can be understood that the dynamics of black hole formation is determined within the background Hubble 
length and is not affected by a global mass excess. Both the amplitude measures S and Cmax are quasilocal quantities 
and do not care about perturbations on scales much longer than the horizon scale at horizon entry. 

On the other hand, we can see that '00,c shows a very different behavior from Sc and Cmax.c- In Table U and Fig. [2l 
we can see that as a increases and hence the transition is smoother, tpQ c increases. The same behavior is also seen 
in Table El and Fig.J3 As a increases and hence the transition is sharper, 0o_c decreases. Based on this behavior, 
Shibata and Sasaki (9| conclude that if the overdensity is surrounded by a low density region, it efficiently collapses. 
We have confirmed that their conclusion is correct if one directly links the small 0o,c to the efficient production of 
black holes, which should applies if the statistical or probability distribution of 0o is regarded as fixed and is centered 
at its unperturbed value. 

To elucidate why neither Sc nor Cmax but 0o,c significantly depends on the behavior of the surrounding region, it 
is helpful to recall that 4' corresponds to the Newtonian potential for the density perturbation in the CMC slicing, 
which can be most clearly seen in Eq. (14.171) with k = 0. The Newtonian potential contains the information of the 
surrounding region in contrast to the averaged density and the compaction function. For example, we can determine 
the averaged density and the compaction function at some radius r from the density distribution inside the sphere of 
radius r. On the contrary, to determine the Newtonian potential and its central value, we need the density distribution 
not only inside but also outside the sphere. Since the PBH formation threshold should be determined by dynamics 
inside the Hubble length, it would be most efficiently described by quasilocal quantities, such as Sc and Cmax, rather 
than 00,c- This is consistent with the suggestion by Young et al. (l8l| . 

At this stage it would be worthwhile to make several comments about the relationship between the arguments 
above and Nakama et al. [T^ . This paper investigated the PBH formation condition for a much wider class of 
initial curvature profiles K{r), which is described by a function with as many as five parameters, basically using the 
method of Polnarev and Musco m- This function mathematically includes (17.31) and also includes profiles which 
are very close to a top-hat, as well as profiles which are gentler than a Gaussian. In addition, that function was 
claimed to include the profiles investigated by Shibata and Sasaki, characterized by Eq. (EH, since the function 
turned out to fit the profiles of Shibata and Sasaki, after being translated into K(r), fairly well with appropriate 
parameter choice. Actually that function was introduced partially in order to express gentler profiles as the ones 
investigated by Shibata and Sasaki, which are realized when a is large. However, this statement about the inclusion 
may seem mathematically inaccurate since [l5l| also restricts attention to exactly compensated profiles, while Shibata 
and Sasaki’s profiles are strictly speaking overcompensated as has been pointed out above. Still, that statement about 
the inclusion is physically justified since whether the initial perturbation is compensated or not does not affect the 
formation of PBHs that much, as has also been pointed out above. 

Though the definition and physical meaning of S, along with the definition of rg as the radius of the overdensity, is 
clear and hence it is convenient, when one tries to discuss the PBH formation condition for more general profiles more 
precisely, it is also useful to introduce another phenomenological parameter characterizing the amplitude of initial 
perturbations instead of S, since S is too sensitive to the information around rg, which is expected not to affect the 
formation of PBHs that much. In [Ts^ a parameter I was phenomenologically introduced, which is similar to S, as 
well as A, which can be interpreted to measure pressure gradient force. It turned out that only these two parameters. 
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characterizing perturbation profiles, are sufficient to describe the formation of PBHs quite well for the generalized 
class of curvature profiles specified by five parameters. In addition, the tendency mentioned above that 5c is larger 
when the transition is sharper due to larger pressure gradient force, is also manifest in terms of these two crucial 
parameters I and A for the generalized profiles. 


VIII. ANALYTIC RESULTS WITH SIMPLIFIED MODELS 

A. Compensated top-hat density model: Physical interpretation 

Here, we will see why '00,c behaves in the opposite way to 5c and Cmax.c- Since our purpose here is to understand 
this behavior qualitatively, we mimic the density profile model (ED by the following simple function: 


fcMC = Cs [0(ro -r) - a ^0(crro - r)] 

where Q{x) is Heaviside’s step function. From Eq. (16.281) . we find 

( (7^(1-cr"3) (0<r<ro) 

K{r)=l - rl)]/r'^ {tq < r < aro) , 

[ 0 (r > OTo) 

where we have put 


( 8 . 1 ) 


( 8 . 2 ) 


/ 4 

■“ 9r2 2-^ 
'o 


Cs. 


(8.3) 


We can see this model is exactly compensated and that rg and ri coincide with each other. For this model, we can 
find 




= ^cylil - a y, 


5 = 


3r 


:Cyyi-a-y. 


3r + 2'^° ov 

On the other hand, by invoking the linear-order approximation in Eqs. (14.3811 and (16.3L we obtain tpQ as 

V’o - 1 = z^'s^ly - o-"^). 


(8.4) 

(8.5) 

( 8 . 6 ) 


We should note that '0o ~ 1 depends on tr in quite a different manner from Cmax and 5. 

Although this model is very simple, it is still difficult to analytically obtain the black hole threshold. The density 
gradient is initially infinite but it becomes finite immediately after the time evolution sets in, so that we have to take 
the balance between the gravitational force and the pressure gradient force into account in this highly dynamical 
system. Here, we make this model more phenomenological noting the intriguing dependence on a. Eor this purpose, 
we introduce three positive constants of order unity, Ci, C 2 and Cg, to simply parametrize the profile dependence as 
well as nonlinearity. 


ci^Csrl{l-a 1), 

(8.7) 

C2-;jCyQ{l — a ^), 

(8.8) 


(8.9) 


These three parameters are unity for the top-hat density model (18.11) . Erom the above equations, we can derive the 
following phenomenological relations: 


; _ci3F-|-21-crW _C 2 3F-|-2~ 

_o^^max — -1 ^max — 0 . 


4 C 2 1 — cr 


C3 8F 1 - cr-3 


C3 6F 


( 8 . 10 ) 
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Here we just assume that the a dependence of 5c and Cniax,c is very weak. Then, is a monotonically increasing 
function of a for cr > 1. The above simplistic analysis explains the qualitative feature of the numerical results shown 
in Table U and Fig. [2] fairly well. In fact, if we choose 

— ~ 2.8, — ~ 3 (8.11) 

C2 C3 

the relations (18.101) agree with the numerical results in Fig.[2]qualitatively, although the numerical results are obtained 
by the model Eq. (17.81) and the Einstein equation is solved fully nonlinearly. The above discussion suggests that the 
formation criterion can be well described by the quasilocal quantities such as Cmax and 5. To translate this criterion 
into the curvature fluctuation ipQ, we need to know the perturbation profile not only in the overdense region but 
also in the surrounding underdense region which may spread to large distances. We should also note that we do not 
expect that the profile-dependence parameters ci, C 2 and C 3 strongly depend on the equation of state because these 
parameters are determined only by the initial density profile. 


B. Uncompensated top-hat density model: Environmental effect 


The physical argument about '0o,c not only explains why neither 6c nor Cmax but '0o,c is significantly affected by 
the behavior of the surrounding region but also suggests that ■i/'o.c is affected by a significant environmental effect 
if the PBH formation results from a perturbation on top of a perturbation of longer wavelength. Such a situation 
is expressed in the top-hat density model by replacing the compensating underdense region with a perturbed region 
with the density perturbation Si . To describe the overdense region of scale rp with the density contrast Ss surrounded 
by the perturbed region of scale aro with the contrast Si, we generalize the top-hat density model given by Eq. (EH) 
as follows: 


2 — ^ 

ScMC = Cs[Q{ro - r) + qQ{aro - r)] , ( 8 . 12 ) 

where q := Si/Ss- K{r) is then given by 

f C'g{l + q) (0<r<ro) 

K{r)=< Cg[{l + q)rl + q{r^-r^)]/r^ {tq < r < aro) , (8.13) 

[ C^s{l + qa^)rl/r^ {r > aro) 

and Cg is given by Eq. (18.31) . So this model has a nonvanishing mass excess and hence is uncompensated except 
for q = —a~^. We can identify tq and aro with the short wavelength s and long wavelength I, respectively, so that 
a = Ijs. The phenomenological relations (I8.10|) should be replaced with 


1 ^cil + qa^^ 
yO ^ ^ , ^max 

4 C2 1 + g 


Cl 3r -I- 2 1 -I- qa'^ ~ 
C3 8 F 1 -I- g 


C2 3r -i- 2 ~ 


(8.14) 


where we have introduced ci, C 2 and C 3 to simply parametrize the profile dependence and nonlinearity. For the exact 
top-hat density model given by Eq. (I8.12p . we identify all of ci, C 2 and C 3 with unity. Note that Si > —1 if we assume 
the density field is everywhere nonnegative. 

We should here note that to define S we invoke averaging over the overdense region. However, if 5; > 0 in this 
model, the radius of the overdense region is not cq but crro. This is a subtle issue in the definition of the averaged 
density perturbation. We here continue to take tq. In fact, if a is much larger than unity, this ambiguity is rather 
in the choice of the averaging length. If we want to discuss the formation of black hole in the scale of ro, we should 
make averaging with the scale of tq . 

Suppose that g > — 1 and that Cmax.c and Sc very weakly depend on the perturbations of longer wavelength. 
Equation (18.141) implies that 'tpo,c is increased (decreased) when the overdense region of short wavelength scale ro 
is surrounded by an overdense (underdense) region of long-wavelength scale aro with cr > 1 and Si > —Sg- The 
larger the absolute value of the density perturbation ratio |g| or the larger the scale ratio of the perturbations a, the 
larger the environmental effect becomes. On the other hand, the requirement K{r)r‘^ < 1 at r = aro implies that 
qa^ (1 + (l)^c^ for the threshold model in the case of tr ^ 1 and hence the environmental effect in ^o.c remains of 
order unity if g > 0. On the other hand, if g < 0, that is, if the overdense region is located in a larger underdense 
region, there seems no limit on the environmental effect, although this must be interpreted with caution because the 
above estimate of ipo only relies on the linear-order analysis. Hence, if one links the small value of ^o,c to the efficient 
production of PBHs, one may conclude that the PBH production is significantly enhanced. It should be noted that 
the above analysis is simplistic and further analysis is necessary in this context. 
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C. Top-hat curvature model: Sharpest transition 


We can get an analytic expression for the top-hat curvature model, where K (r) is given by 

K{r) = 0(A — r), 

with 0 < A < 1. For 0 < r < A, we can find 

^ d- + 

Jo ^ \ a/I - K{x)x^ j 2 

Then, Eqs. (16.211) yield 

1 + y/l-A^ 2(1-f Vl - A2 )cc7 

w = r - , or r = - , - 

1 -I- \J\ — (1 -I- Vl — A^)^ -I- 

for 0 < r < A or 0 < -07 < A, while r = -ccj for r > A or ccj > A, and 

. 1 + VT^ . 2(1 + V1-A2) 

1 + Vl — A^ (1 + Vl — A^)^ + 

for 0 < r < A or 0 < -07 < A, while tj/ = 1 for r > A or ro > A. 

Equations (I4.68|) . (16.151) and (16.351) give Sc, Sc and Ccmc, where Eq. (16.201) gives 

/ = 0(A-r)-y(5(A-r). 


(8.15) 


(8.16) 


(8.17) 


(8.18) 


(8.19) 


The top-hat curvature model is clearly distinct from the top-hat density model because of the delta function in 
the density perturbation. This model is unphysical because the density perturbation S is infinitely negative at the 
transition r = A from the overdense region to the FRW region. However, if the continuous model has a very sharp 
transition, then we expect that the top-hat curvature model may describe the dynamics of the continuous model 
approximately. 

The top-hat curvature model has infinite density gradient and hence infinite pressure gradient force for general 
equations of state. Moreover, since the density perturbation itself is negatively infinite at the transition due to 
the presence of the delta function, the pressure gradient force always dominates gravitational attraction and hence 
prevent the model from collapsing to a black hole except for the case in which the overdense region is initially trapped. 
This implies that the criterion for the black hole formation should be given by K{r)r^ = 1 at the transition, i.e., 
A = 1, because the coordinate singularity at which K{r)r‘^ = 1 implies a marginally trapped surface [T^ . This means 
Sc = 3r/(3r + 2) and Cmax,c = 1/2, respectively. The result of the numerical simulation by Polnarev and Musco [ll| 
for the sharpest top-hat-type profile is consistent with this argument. If we substitute A = 1 and r = w = 0 into 
Eq. (I8.18L we find -tpo = V^. This is close to il)o,c for a not so larger than unity as seen in Table HI The similar trend is 
also seen in Table |TT] for the larger values of a. Since these models have a steeper transition from the overdense region 
to the FRW universe, our expectation is supported by the numerical result. That is, the threshold values for the model 
with a sufficiently steep transition will be approximately given by those for the top-hat model. Sc = 3F/(3F -|- 2), 
llmax,c — 1/2 and ^o,c — \/2- 


D. Double top-hat curvature model: Environmental effect 

Next we move onto a double top-hat curvature model, for which K{r) is given by 

K{r) = (1 - H)0(Ai - r) + H0(A2 - r), 


( 8 . 20 ) 


where A is constant, 0 < Ai < A 2 and we require K(r)r^ < 1 for the avoidance of the coordinate singularity. A 
similar model is used in Nakama [l^ in a different context [s^- We now extract an environmental effect from this 

simple model. For this model, Eqs. (14.681) . (I6.15|) and (16.351) give Sc, Sc and Ccmc, where Eq. (j6.20l) gives 


/ = ( 1 -^) 


0(Ai-r)-^5(Ai-r) 


0 (A 2 - r)- ^i5(A2 - r) 


( 8 . 21 ) 
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and from Eq. (|6.18l) we can find 


"00 


2(1 + \/l — 


(1 + \/r^^Af)(l + \/l — A^A) 


( 8 . 22 ) 


We can see that the single top-hat curvature model is reproduced by putting ^ = 0 or A = 1 or Ai = A 2 . We can 
make identification s = Ai, Z = A 2 and q = 61 /Ss = A/— A), and hence the central overdense region is surrounded 
by the lower but still overdense region if 0 < A < 1, while the surrounding region is underdense if ^ < 0. We should 
however note the delta function contributions in the density profile again. 

Similarly to the single top-hat curvature model, the threshold for black hole formation in this model is given 
by max(A^,A 2 A) = 1. Let us focus on the cases where A|A < 1. Then, Ai = 1 gives a threshold and hence 
(5c = 3r/(3r -I- 2) and Cmax.c = 1/2. Therefore, 5c and Cmax.c do not depend on the surrounding region. However, this 
is not the case for which is calculated as 


00,0 = 


2(1 + VT^) 

1 + ^1-AjA' 


(8.23) 


Figure S] shows 0o,c as a function of A 2 for different values of A. Noting A 2 > Ai = 1 and A^H < 1, we can conclude 
that ipQ^c > for H > 0, while 00,c < for H < 0. If A 2 (> 1) is fixed, 0o,c monotonically increases from y^2/A2 to 

Y^^2(T-i-~y/r^T7^f) as A increases from —00 to 1/A|. If H £ (0,1) and A is fixed, 0o,c monotonically increases from 

\/2 to 2(1 -I- -01 — A) as A 2 increases from 1 to Xj'/A. If H < 0 and A is fixed, tjjo,c monotonically decreases from 

y/2 to 0 as A 2 increases from 1 to 00 . Thus, the existence of the surrounding overdense (underdense) region of longer 
wavelength increases (decreases) 0o,c and, hence, suppresses (enhances) the PBH production if one assumes that the 
smaller 0o,c is directly related to higher production rate. Moreover, the environmental effect on 0o,c is bounded if the 
density perturbation of longer wavelength is positive, while 0o,c even gets smaller than the unperturbed value if the 
wavelength I of the underlying negative density perturbation is sufficiently long and q = Si/Ss is fixed, and approaches 
0 in the limit where I is infinitely long. The qualitative behavior of 0o,c is common for both the uncompensated 
top-hat density model and the double top-hat curvature model, in spite of the physical difference between the models. 



FIG. 4. The threshold values of ipo.c for the double top-hat curvature model given by Eq. (18.2011 as functions of A 2 with different 
values of A. 


IX. CONCLUSION 

We have constructed cosmological nonlinear perturbation solutions with the long-wavelength scheme in the CMC, 
uniform-density, comoving and geodesic slicings without assuming symmetry. These solutions are generated by only a 
master variable fk = where the perturbations from the flat FRW solution can be arbitrarily large. We have also 
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derived the explicit relation among these four slicings. One of the interesting applications is to construct initial data 
for primordial structure formation with and without spherical symmetry. For example, we can study nonspherical 
effects on the PBH formation, which are expected to be important especially for the soft equation of state. 

Then, we have presented two distinct formulations of spherically symmetric spacetimes and the definitions of mass 
excess, compaction function and averaged density perturbation as spatial gauge invariant but slicing dependent per¬ 
turbation quantities. Based on the general formulation of long-wavelength solutions, we have constructed spherically 
symmetric long-wavelength solutions in the CMC slicing and in the comoving slicing. We have elucidated the relation 
between the two spatial coordinates, the conformally flat coordinates and the areal radial coordinates. Using these 
solutions and relations, we have established the equivalence between the long-wavelength solutions given in Shibata 
and Sasaki Q and the asymptotic quasihomogeneous solutions given in Polnarev and Musco [llj | , both of which have 
been used as initial data sets for the simulations of PBH formation. 

Using this equivalence, we have reproduced the numerical simulation by Shibata and Sasaki Q with the numerical 
code based on the formulation by Polnarev and Musco m and obtained the results which agree with the result 
of Shibata and Sasaki Q. We have also calculated ijjo,c for the numerical simulation by Polnarev and Musco m- 
Combining these results, we have discussed that the smoother the transition from the overdense region to the FRW 
universe is, the smaller the 6 c and Cmax.c become. The minimum values of Sc and Cmax.c are attained if the transition 
from the overdense region to the homogeneous universe is the smoothest. We have discussed that the analytic 
threshold formula obtained by Harada et al. [T^ should apply for this case, while the analytic formula for the possible 
maximum value can be regarded as the threshold value for a sufficiently sharp transition. We have found the relation 
among ipo^ Cmax and 6 by using the compensated top-hat density model and interpreted why '0o,c shows an apparently 
opposite behavior within the (approximately) compensated models to those of dc and Cmax.c- Moreover, generalizing 
the top-hat density model to the uncompensated one, we suggest an environmental effect on '0o,c from the underlying 
long-wavelength perturbations. This is also supported by the double top-hat curvature model, which can be treated in 
a fully analytic and nonlinear manner. We conclude that even if 6 c and Cmax.c are not sensitive to the density profiles 
in the surrounding region, i/jq^c can be significantly larger (smaller) if the density perturbation of the wavelength in 
which we are interested is in the underlying positive (negative) density perturbation of much longer wavelength. 
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